Fast transitionless expansions of cold atoms in optical Gaussian beam traps 



E. Torrontegui/ Xi Chen,^''^ M. ModugnOj^^"* A. Ruschhaupt,^ D. Guery-Odelin,^ and J. G. Muga^-"^ 

Departamento de Qmmica Fisica, Universidad del Pais Vasco - Euskal Herriko Unibertsitatea, Apdo. 644, Bilbao, Spain 
'^Department of Physics, Shanghai University, 200444 Shanghai, P. R. China 
'^Departamento de Fisica Teorica e Histona de la Ciencia, 
Universidad del Pais Vasco - Euskal Herriko Unibertsitatea, Apdo. 644; Bilbao, Spain 
^ IKERBASQUE, Basque Foundation for Science 
^Institut fur Theoretische Physik, Leibniz Universitdt Hannover, Appelstrafie 2, 30167 Hannover, Germany 
" Laboratoire Collisions Agregats Reactivite, CNRS UMR 5589, IRSAMC, 
Universite Paul Sabatier, 118 Route de Narbonne, 31062 Toulouse CEDEX 4, France 
''Max Planck Institute for the Physics of Complex Systems, Nothnitzer Str. 38, 01187 Dresden, Cermany 

We study fast expansions of cold atoms in a three-dimensional Gaussian-beam optical trap. Three 
different methods to avoid final motional excitation are compared: inverse engineering using Lewis- 
Riesenfeld invariants, which provides the best overall performance, a bang-bang approach, and a 
fast adiabatic approach. We analyze the excitation effect of anharmonic terms, radial-longitudinal 
coupling, and radial-frequency mismatch. In the inverse engineering approach these perturbations 
can be suppressed or mitigated by increasing the laser beam waist. 
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I. INTRODUCTION: DRIVEN EXPANSIONS OF 
COLD ATOMS 

Driving the expansion of a Bose-Einstein condensate 
or more generally of a cold atom cloud in a controlled 
way, by implementing a designed time-dependence of the 
trap frequency, is a common and basic operation in a 
cold-atom laboratory. The expansion may be aimed at 
different goals, such as decreasing the temperature 0, [2] , 
adjusting the density to avoid three body losses m, fa- 
cilitating temperature and density measurements S^, or 
changing the size of the cloud for further manipulations 
Expansions, isolated [5] or as part of refrigeration 
cycles 16|, are also important from a fundamental point 
of view to quantify the third principle of thermodynam- 
ics.^ In general, changes of the confining trap will excite 
the state of motion of the atoms unless they are done 
very slowly or, in the usual quantum-mechanical sense of 
the word, "adiabatically" , but slow-change processes are 
also prone to perturbations and decoherence, or imprac- 
tical for performing many cycles. Engineering fast ex- 
pansions without final excitation is thus receiving much 
attention recently both theoretical and experimental 
[22j . Most theoretical treatments so far are for idealized 
one-dimensional (ID) systems, but the implementation 
requires a three-dimensional (3D) trap [ll|, [13, in 
principle with anharmonicities and couplings among dif- 
ferent directions. 

In this paper we address these important aspects to im- 
plement the expansion in practice. Specifically we shall 
model a simple physical realization based on an elongated 



Compressions may also play a role in many experiments but their 
treatment is similar to expansions so we shall not deal with them 
here explicitly. 




FIG. 1: (Color online) Schematic view of the optical trap. 

cigar-shaped optical dipole trap with cylindrical symme- 
try, see Fig. [TJ This trap is formed by a single laser 
which is red detuned with respect to an atomic transi- 
tion to make the potential attractive (Section |TT]), and is 
characterized in the harmonic approximation by longitu- 
dinal and radial frequencies. While magnetic traps allow 
for an independent control of longitudinal and radial fre- 
quencies [III, [13, [13 J this is not the case for a simple laser 
trap that therefore requires a special study. We assume 
that the time-dependence of the longitudinal frequency 
is engineered to avoid final excitations with a simple ID 
harmonic theory (Section IIII[) and analyze the final fi- 
delity in the actual trap. Even though for full 3D-results 
we resort to a purely numerical calculation in Section Ivl 
an understanding of the effects involved is achieved first 
by analyzing separately longitudinal and radial motions 
in Section ITVl Conclusions and open questions are drawn 
in the final Section IVTl 

II. THE MODEL 

The intensity profile of a Gaussian laser beam in the 
paraxial approximation is given by 

/(r,z,^) = /o(Oe-^-^/-^(-) \, (1) 

where r and z are the radial and longitudinal coordinates 
respectively, and the variation of the spot size w with z 
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is given by 



w{z) 




where zr — ttwq/X is the Rayleigh range, wq the waist, 
and A the laser wavelength. The paraxial approximation 
is valid for waists larger than about 2A/7r. 

If an atom is placed in this field with a detuning 5 
which is large with respect to the Rabi frequency ft and 
the inverse of the lifetime of the excited state F = t~^, 
but small with respect to the transition frequency, inter- 
nal excitations and counter-rotating terms are negligible, 
and the potential that the ground state atoms feel due 
to the dipole force is proportional to the laser intensity. 
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and inversely proportional to 6. Isat = ttHc/ (SX^t) is the 
saturation intensity.^ Combining Eqs. ([T]) and ([2]), and 
adding for convenience in later expressions the physically 
irrelevant term Vq (t) , the potential takes finally the form 

V{r,z,t) = -Vo{t)e-''''^'"'^' \ \ , +V,{t), (3) 

i + Z I Zj^ 

where Vo(t) = lQ{t)hT'^ / {SSIsat}- The analysis carried 
out here is not restricted to a two-level atom. Indeed, 
the dipole force is always proportional to the intensity 
for a multi-level atom but with a coefficient that depends 
on the level structure and the light polarization [23] • As 
an example of the magnitudes involved in practice, for 
rubidium-87 atoms, a Gaussian beam linearly polarized 
of waist 8 /im, power P = 30 W and wavelength 1060 nm 
provides a longitudinal trap frequency uj^^jli: ~ 2500 
Hz. We shall use these or similar values in numerical 
examples below. 

In this work we shall assume small enough densities 
and times so that a one-particle, Schrodinger descrip- 
tion for the atomic motion neglecting collisions is valid. 
Extensions to Tonks-Girardeau gases or condensates are 
possible as in 0, H, . We shall also assume that the 
longitudinal axis lies horizontally and neglect the effect 
of gravity, either because it is artificially compensated, 
or because the radial confinement is tight. 



^ Indeed other regimes could be used as well 23] . For larger detun- 
ings counter-rotating terms become important, but the dipole po- 
tential is still proportional to the laser intensity. In the extreme 
limit of quasi-electrostatic traps the frequency of the trapping 
light is much smaller than the resonance frequency precluding 
the possibility of a potential sign change from attractive to re- 
pulsive; this change is in principle allowed by Eq. l(2]| by changing 
the sign of the detuning across resonance, but see the discussion 
below. 

^ This requires that the vertical shift of the trap center in the 
gravity field is small with respect to the waist, or g <^ wgui-^, 
where loh is the radial angular frequency, and g the acceleration 
due of gravity. 



To solve the time dependent Schrodinger equation as- 
sociated with the potential in Eq. ([3]), 



at 2m 



(4) 



we use cylindrical coordinates. 
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With the standard separation of variables ^ = 
J^{r,z,t)^{(f>), we find 
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(the simple or double subscripts denote first or second 
derivatives) whereas the equation for $ reads 



2 ^ rj. 



(5) 



where v is a, magnetic quantum number determining the 
conserved angular momentum component along the z di- 
rection. This is solved by e*"*^ so we shall concentrate on 
J-{r, z, t). Defining now ^'(r, z, t) = \/rJ-{r, z, t) we get 
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where 
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+ y(r,z,i). 



(6) 
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and V{r,z,t) is given by Eq. ([3|). The task is now to 
design Vo{t) so as to achieve a fast expansion without 
final excitation. 



III. EXPANSION PROTOCOLS 

In this section we briefly present three expansion pro- 
tocols for a one-dimensional harmonic trap. They are 
chosen mostly because they have been realized experi- 
mentally [H, [nl, [13, [13] • Their simplicity will help us 
to identify relevant physical phenomena that will affect 
also other protocols. No claim of completeness or global 
optimization is made, and there is certainly room for in- 
vestigating other protocols. 

The objective is to expand the trap from an initial fre- 
quency a;o/27r to a final frequency a;//27r in a time tf, 
without inducing any final excitation in the adiabatic or 
instantaneous basis for which the Hamiltonian is diago- 
nal. Note that transient excitations are generally allowed. 
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A. Inverse engineering with Lewis-Riesenfeld 
invariants 

Lewis-Riesenfeld invariants may be used to engineer 
efficient expansion protocols as explained in 0, H, [HI . 
We refer the reader to these works for the details and 
give here only the elements required for a practical ap- 
plication. For any harmonic oscillator expansion, and in 
fact for any potential with the structure 



(SI) 



where the function U is arbitrary, there is a quadratic- 
in-momentum invariant operator of the form 



I(q,p,t) = - — [b^p'^ — mbb{qp + pq) + m^b'^q^] 
Zm 



1 2 /9 
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(10) 



where q and p are generic position and momentum op- 
erators (particularized later for longitudinal or radial di- 
rections), and b = b(t) is a scaling function that satisfies 
the Ermakov equation 



b + uj\t)b^uj^/b-^. 



(11) 



To inverse engineer the trap frequency, b{t) is designed to 
satisfy the boundary conditions that guarantee no final 
excitations and continuity of the trap frequency. 



6(0) = 1, 6(0) = 0, 6(0) = 0, 
b{tf) = 7, bitf) = 0, bitf) = 0, 



(12) 



where 7 = y^ojQ/ujf. Finally w(t) is deduced from Eq. 
([TT]) . When these conditions are satisfied the n-th eigen- 
state of H{0) ^ p'^/{2m) + V{q,0) evolves as an "ex- 
panding mode" , which is the time-dependent n-th eigen- 
vector of the invariant times a "Lewis-Riesenfeld" phase 
factor e'"", where a„(i) ~ —{n + l/2)ajo /J dt' /b^, until 
it becomes an n-th eigenvector of the final Hamiltonian 
H{tf)=py{2m) + V{q,tj). 

Discontinuities in 6 may in principle be allowed, but 
they amount to perform finite jumps of the trap fre- 
quency. This is an idealization but it can be approached 
up to technical limits. Discontinuities in 6 have more se- 
rious consequences, as they imply delta functions for 6 
and infinite trap frequencies at the jumps. Thus their 
approximate physical realization is even more difficult, 
but they are useful to find bounds for physical variables 
of interest as in [5| or in Appendix [X] There are many 
b{t) that satisfy Eq. ([T^. A simple smooth function is a 
quintic polynomial with six coefficients, 

b{t) =6(7- l)s''^-15(7- 1)/ + 10(7- 1)/ + 1, (13) 

where s := t/tf. More sophisticated choices are possible 
to minimize the process time, the average energy, or other 



variables with or without imposed constraints ( "bounded 
control" ) on the allowed trap frequencies [H, [11] . 

In principle this method allows for arbitrarily small val- 
ues oi tf but for very small process times the potential 
becomes transitorily repulsive ^S]. With a laser beam, a 
repulsive potential could be achieved by means of a pos- 
itive (blue) detuning. However the validity of Eq. Q 
requires a large detuning, so care should be exercised 
when crossing the resonance. In this work the numerical 
examples are restricted to positive trap frequencies since 
the transient passage from attractive to a repulsive inter- 
action could involve unwanted effects such as radiation 
pressure. We shall in any case point out in the figures 
the minimum value oi tf for which the potential stays 
attractive for all times with the quintic b{t). 



B. Bang bang 

Bang-bang methods are those with a stepwise constant 
behavior of some variable. The simplest one to avoid final 
excitations consists of using a constant intermediate trap 
frequency, the geometric average of the initial and final 
frequencies, 

f wo, t<0 

Jang ^ J (^0,^^)1/2^ Q<t< t^""^ 

(14) 

during a fourth of the corresponding period. Optimal 
bang-bang trajectories with two intermediate frequencies 
and arbitrary final times are also possible 0, [^ but they 
are not considered here. To arrive at Eq. (jl4p . one may 
use a classical argument [13] or proceed as follows: the 
solution of the Ermakov equation for t > assuming a 
constant intermediate angular frequency wi with initial 
boundary conditions 6(0) = 1, 6(0) = is 

b{t) = ^[{loI-loD/loI] sm\uj,t) + l. (15) 

Then, if we solve for t f and ui the two equations implied 
by the final boundary conditions 

6(i;)=7, 6(t/)-0, (16) 

the values in Eq. are found. 

C. Fast adiabatic protocols 

"Fast" and "adiabatic" may appear as contradictory 
concepts. A fast adiabatic protocol seeks to perform a 
process as quickly as possible keeping it adiabatic at all 
times. For harmonic oscillator expansions the adiabatic- 
ity condition reads w/w^ ^ 1, so making wjuP' constant 
from < = to i/ , l| , and solving the resulting differential 
equation for a;(i), we get [8| 



;°*(t) = 



l-(c^y-wo)^' 



(17) 
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IV. LONGITUDINAL AND RADIAL MOTIONS 



By expanding the potential in Eq. 
Taylor series around (z = 0, r ~ 0), 



([3]) in a double 
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we see that the first coupling term between radial and 
longitudinal motions is of fourth order, proportional to 
r^z^. If we could neglect this and higher order coupling 
terms, longitudinal and radial motions would be indepen- 
dent. The approximation that considers longitudinal and 
radial motions completely uncoupled is useful to analyze 
several effects separately and gain insight. Not only that, 
as we shall confirm later with full 3D calculations, this 
is also a good approximation numerically for low energy 
levels. 



A. Longitudinal motion 

To study the motion in the longitudinal direction we 
consider first the full longitudinal Hamiltonian (putting 
r = in Eq. ©) 



H{z,t) 



2m dz^ 



Vo{t) 



1 



1 + zV4, 



(19) 



To characterize this potential in terms of a longitudinal 
frequency consider the harmonic approximation 



Hhar{z,t) — — 



2m dz^ 



52 , , z2 

+ Vo{t)- 



Pz 

2m 



mLo1{t)z'^ 



where Pz = —ihd/dz and 



2Vb(0 



(20) 



(21) 



We may now apply the methods described in Sec. IIIII 
We shall add a subscript z, for "longitudinal", to the 
trap frequencies and take q ^ z. For an imposed uJz{t) 
and fixed waist and laser frequency, Eq. (I?T|) determines 
the time dependence of the laser intensity. 

Figure [5] shows for the ground state, n = 0, the "lon- 
gitudinal fidelity" Fl - |(Z„(t/)|C/^(t/,0)|Z„(0))|, where 
Zn{z, 0) and Zn{z,tf) arc the initial and final n-th eigen- 
states of the full longitudinal Hamiltonian, Eq. ([TS]). 
for frequencies ujoz and tufz respectively, and Uz{tf,0) 
is the evolution operator with Eq. (IT9)) . As in all figures 
hereafter the initial longitudinal frequency of the trap 
is Ldoz/2Tr — 2500 Hz, and the wavelength is taken as 
A = 1060 nm, characteristic of Neodymium-doped lasers. 
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FIG. 2: (Color online) Longitudinal fidelity for Eq. (HH) ver- 
sus final time tf for two different final frequencies: in (a) 
iUfz/2n = 250 Hz; in (b) ijjfz/2n — 25 Hz. In both figures 
(jjoz/'2tv — 2500 Hz, A = 1060 nm, and the initial state is the 
ground state. All curves have been calculated for two diflferent 
waists, Wo = 10 fim and wq — 3 ^im, but the results are indis- 
tinguishable. Vo{t) is chosen according to the three protocols: 
inverse engineering, using Eq. (fT3)) in Eq. (fTTjl (solid red line); 
bang-bang (circles for wq = 10 fim and triangles for uiq = 3 
fj,m); and fast adiabatic method (long-dashed blue line). The 
green vertical lines denote the minimum tf for which Vo{t) 
remains positive for all t in the inverse engineering approach 
with a quintic b{t). 



The vertical green lines in this and the following figures 
correspond to the minimum tf value for which ujz{t) is 
positive for all t using the quintic b{t). Of course this 
limit only applies to the inverse engineering approach. 

The three methods are compared for two different fi- 
nal frequencies, in (a) ujfz/2Tr — 250 Hz and in (b) 
ujfz/2Tr — 25 Hz, and two different waists, see the caption 
for details. Actually the effect of the waist change is neg- 
ligible in the scale of these figures. Globally the inverse 
engineering method outperforms the others. The bang- 
bang approach provides a good fidelity, but only for a 
specific final time, and the fast adiabatic method fails in 
fact to be adiabatic at short times. This is very evident 
for the smaller final frequency in Fig. [2] (b): the opening 
of the trap is faster and the level spacings smaller than 
for the larger final frequency, so the state cannot follow 
an adiabatic behavior. 
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FIG. 3: (Color online) Longitudinal fidelity versus level num- 
ber n. For Wo = 3 jj,m (blue symbols): fidelity Fl ~ 
\{Zn{tf)\Uz{tf,0)\Zn{0))\ with inverse engineering computed 
for Eq. I|19|l and a quintic b{t) (empty triangles); first-order 
bound in Eq. (|A9|) (x); second-order approximate fidelity 
computed with Eqs. (|A11[) and (jA12l) for the quintic b{t) 
(filled triangles). For a larger waist, wq — 10 fim (red sym- 
bols), the corresponding values (empty circles, crosses, and 
filled circles), are very nearly one and indistinguishable. In 
all cases a;/z/27r — 25 Hz, ljoz/2-k = 2500 Hz, A = 1060 nm, 
and tf = 2.5 ms. 



For small tf and/or large 7 the transient excitation 
energy during the inverse engineering protocol may be 
high, giving the atoms access to the anharmonic part of 
the potential. This could lead to a decay of fidelity as tf 
decreases, but the effect is not seen in the scale of Fig. 
[21 The small effect of anharmonicity is enhanced by in- 
creasing the vibrational number, see Fig. [3] It can be 
avoided by increasing the waist, which reduces the anhar- 
monic terms, as demonstrated in Fig. [3] The Appendix 
[^provides a perturbation theory analysis of longitudinal 
anharmonicity. A first order approach gives analytical 
lower bounds for the fidelity, and the possibility to maxi- 
mize them by other choices of b{t). More accurate results 
are found in second order, see Fig. [3l 



B. Radial motion 

To study radial motion we define the radial Hamilto- 
nian by setting z = in Eq. ([3]) , and add the "centrifugal 
term" (an attractive one for v — 0), 



H{r,t) = --^^^l/o(<)(e-2'-'/-o_i 



ff_ 1^1/2 _ 1/4 
2m 



In the harmonic approximation we take 

muj'j^{t)r'^ 



(22) 



(23) 



_ 4Fo(i) 



(24) 



defines the radial frequency of the trap uj nit) / (2Tr) . 

We assume that Vq (t) is set to satisfy the designed lon- 
gitudinal expansion according to Ec^. (|2T|) . Substituting 
this into Eq. ([M)) gives the relation between radial and 
longitudinal frequencies. 



V2: 
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(25) 



which is key to understand the behavior of the radial 
wavefunction. The waist value wq = X/{y/2n) would 
make both frequencies equal but for such a small waist 
the paraxial approximation fails and the present theory 
could not be applied [26j . 

Now we define the ground state "radial fidelity" as 
Fr = \{Mtf)\UR{tf,0)\cbomi where [//?(*/, 0) is the 
radial evolution operator for H{r,t), and 0o(r, 0) and 
(j)o{r,tf) are ground states of the initial and final radial 
traps for Eq. ([22]), with Vo{t) given by Eq. In 
Fig. m the radial fidelity is depicted for the three pro- 
tocols explained above and v = 0. The bang-bang tra- 
jectory causes much excitation and gives a rather poor 
fidelity in the radial direction (see the big empty sym- 
bols). This is because, even though the "correct" tran- 
sient radial frequency a;^°"^(i) = {\/2ttwo / X)uj'^°'"^ (t) is 
implemented (as the geometric average of initial and fi- 
nal radial frequencies), the time i/ is adjusted for the 
longitudinal, not for the radial frequency, compare t^™^ 

to ty^<' := Xt%"y{V2TTWo). 

The behavior of the other two methods is much more 
robust as shown by their fidelities in Fig. 01 where the 
inset amplifies the smal\-tf region. The fast adiabatic 
approach, blue long-dashed lines of Fig. 01 gives an ex- 
cellent radial fidelity. A detailed calculation shows that 
the adiabaticity condition for the radial Hamiltonian, 
see the Appendix |B1 takes, in spite of the centrifugal 
term, the same form as for the ordinary harmonic os- 
cillator, namely ujR/ujj^ <C 1. If lOz/loI is constant, the 
ratio ujr/uj'^ will be constant too, but smaller by a factor 
\/{\/27r'Wo), so that the radial dynamics will be in gen- 
eral more adiabatic. Indeed, comparing blue long-dashed 
lines of Figs. [5] and |H we see that the radial fidelity is 
better than the longitudinal one, as the confinement is 
tighter radially than longitudinally, and the energy lev- 
els are more separated. 

For the inverse engineering protocol with a quintic 6, 
red solid and black short-dashed lines in Fig. 21 the radial 
excitation is small, the fidelity being nearly one for tf > I 
ms in the scale of the figure. This could be surprising 
since the radial frequency does not behave according to 
an ideal inverse engineering dependence 
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FIG. 4: (Color online) Radial fidelity for different final times. 
Parameters: 1^ = 0, LJoz/2n = 2500 Hz, ujf2/2n = 250 Hz, 
and A — 1060 nm. Protocols: inverse engineering, quintic b, 
Wo ~ 10 fim (solid red line); inverse engineering, quintic b, 
Wo = 3 fim (short-dashed black line); bang bang, wq = 10 fim 
(empty circle); bang bang, wo = 3 /im (empty triangle); fast 
adiabatic, wo = 10 pm (long-dashed blue line); fast adiabatic, 
Too ~ 3 fj.m (dotted orange line). In the inset, corresponding 
(red) triangles, (black) squares, (blue) circles, and (orange) 
diamonds are the approximate fidelities from adiabatic per- 
turbation theory using Eq. (I28p . The green vertical line near 
0.4 ms is the minimal final time tf for which the quintic-& in- 
verse engineering protocol implies positive frequencies at all 
transient times. 



where we have used the Ermakov equation and 7 — 
(woz/w/z)^/^ = (ajoi?/w/_R)^/^, so the function b is the 
same for the longitudinal and radial directions. Instead, 
the actual radial frequency varies according to Eq. (1^ , 

The difference between Eqs. ([26l) and ([27l) may be quite 
significant, as the example in Fig. [5] shows, so the ra- 
dial state does not really follow the expanding modes 
corresponding to Eq. ([26l) . As a consequence a pertur- 
bative approach based on Eq. (j26p as the zeroth order 
fails completely, even for qualitative guidance. (Further 
evidence is provided in Fig. [S]) Instead, a valid pertur- 
bation theory analysis may be based on the zeroth order 
of the adiabatic modes, which are instantaneous eigen- 
states of the radial Hamiltonian. The physical reason for 
the relatively good radial behavior of the invariant-based 
inverse engineering protocol is thus the adiabaticity in 
that direction. When tf approaches the critical lower 
value, ujR{t) becomes too steep at small values of w^j, see 
Fig. [SJ and adiabaticity breaks down, which explains the 
fidelity decrease there. To mitigate this problem, a larger 
waist may be used to increase ujR{t) for a given uJzit), see 
Eq. ([25)) . making the radial process more adiabatic and 
improving the fidelity, as shown in Fig. 21 Keep in mind 
that to implement a given uJz (t) an increase of waist has 
to be compensated by an increase of laser intensity. 
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FIG. 5: (Color online) (Square of) radial frequency versus 
time: Eq. (|27|) based on a quintic b (red solid line); fast 
adiabatic protocol, Eq. (I25p with Eq. (|17p for uj^ (blue long 
dashed line); Eq. (|26|) (black short dashed line). Parameters: 
A = 1060 nm, Luoz/2n = 2500 Hz, LUfj2n = 250 Hz, toq = 3 
fj,m and tf = 0.36 ms. 

Inserting a wave function expansion in an adi- 
abatic basis (instantaneous eigenstates) ip{r,t) — 
Y.kak{t){r\(j)k{t))e-'^oEu(t')dt' ^ -j^^Q ^l^g 

Schrodinger equation provides a set of coupled dif- 
ferential equations. Integrating formally for Ofe(t) and 
assuming that in zeroth order af\t) = So,k, we get in 
first order [13,111 

a^^\t) - - rdi'(0i(t')|0o(t'))e-^^°'''*"'''''^*"^~''^'*"'l 
Jo 

= _ rrfi'ii^e2^/o''i*"-«(*"), (28) 
Jo 2wfl(t') 

where the second line follows by applying, in addition, 
the harmonic approximation. All other first order am- 
plitudes (for k > 2) are zero so the radial fidelity at 
tf can be calculated in this approximation as Fji{tf) = 

[1— \a^i\tf)\'^]^^'^. This is in fact correct up to second or- 
der and it provides good agreement with the exact results 
as shown by the symbols in the inset of Fig. |4l 

V. FULL 3D ANALYSIS 

Finally we study numerically the actual, coupled dy- 
namics driven by the exact full potential, Eq. ([3]), 
combining all the effects considered so far and the 
longitudinal-radial coupling. The 3D fidelities are also 
compared with simpler ID fidelities. 

Before discussing the results a remark on the technical- 
ities of the numerical calculations of the time-dependent 
wave functions is in order. To solve the longitudinal time 
dependent Schrodinger equation we approximate the evo- 
lution operator Ul using the Split Operator Method 
(SOM) f29( and Fast Fourier Transform (FFT) 
For the radial time dependent Schrodinger equation, we 
have used instead the finite differences Crank-Nicholson 
scheme [sH. In this way the singular point r = can 
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FIG. 6: (Color online) (Modulus of) the overlap integral be- 
tween the state evolving from the ground state with the radial 
Hamiltonian, using the actual radial frequency in Eq. (|27|l . 
and (a) the instantaneous eigenstate in harmonic approxima- 
tion (red solid line), or (b) the expanding mode in harmonic 
approximation for the radial frequency in Eq. (|26[) (dashed 
blue line). Parameters: u — 0, uoz/Stt — 2500 Hz, A = 1060 
nm, ijjfz/2n = 250 Hz, wq — 3 /im, t/ = 0.6 ms. 
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FIG. 7: (Color online) Fidelity for the 3D ground state versus 
different final times. Parameters: i/ — 0, ijjoz/2tt = 2500 
Hz, A = 1060 nm, ujfz/^n — 250 Hz. Protocols: inverse 
engineering, quintic b, wo = 10 ^im (red solid line) , and radial 
fidelity (red filled triangles); inverse engineering, quintic b, 
Wo = S fim (black short-dashed line) , and radial fidelity (filled 
squares); bang bang, wo = 10 /im (empty circle in inset), and 
radial fidelity (filled circle); bang bang, u;o = 3 jj,m (empty 
triangle in inset), and radial fidelity (filled triangle in inset); 
fast adiabatic, liio = 3 fim and wo = 10 (blue long-dashed 
line), and longitudinal fidelity (blue diamonds). The green 
vertical line marks the threshold for positive trap frequencies 
in the inverse engineering, quintic-fe protocol. 



be excluded imposing a zero of the wavefunction there 
for all t. For the 3D calculation we use similarly the 
Crank-Nicholson scheme with operator splitting for mul- 
tidimensions (the radial and longitudinal directions) [3l| . 

We calculate, starting from the ground state of the 
initial trap for j/ = 0, the fidelity of the final state with 



respect to the ground state of the final trap using the 
three expansion protocols. Figure [7] clearly shows that 
the inverse engineering protocol provides the best over- 
all performance, with a fidelity decay at smaller times 
that can be avoided by increasing the waist. The radial- 
longitudinal quartic coupling does not play any signifi- 
cant role for the parameters of the example, and in any 
case it may also be suppressed by a waist increase. 

As for the simple bang- bang approach, it fails in 3D, 
as expected, due to its poor radial behavior. (The radial 
fidelity alone is close to the 3D fidelity.) The results for 
fast adiabatic and inverse engineering methods are not 
too different from each other for the chosen expansion, 
but the fast adiabatic method will fail sooner for more 
demanding expansions with smaller final frequencies, as 
clearly illustrated in Fig. [l](b), due to its inability to re- 
main really adiabatic for the longitudinal direction. The 
main limitation of the inverse engineering method is in- 
stead the possible failure of adiabaticity in the tighter 
radial direction, so it is intrinsically more robust than 
the fast adiabatic one. Note that the 3D fidelity is mim- 
icked accurately by the ID radial fidelity for inverse en- 
gineering and by the ID longitudinal fidelity for the fast 
adiabatic approach. 



VI. DISCUSSION AND OUTLOOK 

In this work we have analyzed the practical implemen- 
tation of a transitionless expansion in a simple Gaussian 
optical dipole trap, taking into account anharmonicities, 
radial-longitudinal couplings, and the radial-longitudinal 
frequency mismatch. The main conclusion of the study 
is that the transitionless expansions in optical traps are 
feasible under realistic conditions. Despite the relation 
between the longitudinal and transversal trapping fre- 
quencies through the intensity, the different timescales 
enable us to design fast expansions with high fidelities 
with respect to the ideal results using the invariant-based 
inverse engineering method, which is particularly suitable 
compared to the two other approaches examined. Our 
detailed analysis of radial and longitudinal motions re- 
veals the weakest points of each approach: for the inverse 
engineering, the main perturbation is due to the possi- 
ble adiabaticity failure in the radial direction, which can 
be suppressed or mitigated by increasing the laser waist. 
This waist increase would also reduce smaller perturb- 
ing effects due to longitudinal anharmonicity or radial- 
longitudinal coupling. The simple bang-bang approach 
fails because the time for the radial expansion is badly 
mismatched with respect to the ideal time, and the fast 
adiabatic method fails at short times because of the adi- 
abaticity failure in the longitudinal direction. 

Complications such as perturbations due to different 
noise types, and consideration of condensates, gravity ef- 
fects, or the transient realization of imaginary trap fre- 
quencies are left for separate works. Other extensions 
of the present analysis could involve the addition of a 
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second laser for further control of the potential shape, 
or alternative traps. Optical traps based on Bessel laser 
beams, for example, may be interesting to decouple longi- 
tudinal and radial motions. Transport protocols are also 
amenable to a similar 3D analysis. Finally, the combina- 
tion of invariant-based inverse engineering with optimal 
control theory to optimize the frequency design is also a 
promising venue of research 12, 25]. 
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Appendix A: Perturbative treatment of longitudinal 
anharmonicity 

To quantify the effect of anharmonicity in the inverse 
engineering method let us first expand the longitudinal 
Hamiltonian retaining the quartic term. 



V{z) 



1 + z-^/zj^ 



-mujjt)z 



Vo{t) « Vo- 



zi 



:=Vi 



Vo- 



(Al) 



The harmonic term is the unperturbed potential and the 
quartic term the perturbation. A further simplification 
to evaluate the fidelities Fl = |(Z„(t/)|J7^(t/,0)|Z„(0))| 
is to substitute the initial and final exact eigenstates 
l-^n(O)) and \Zn{tf)) by corresponding eigenstates of the 
harmonic oscillator \ipn{0)) and \ipnitf)), see the expres- 
sions below. This is an excellent approximation for the 
lower eigenstates. In Fig. [3l for example, the over- 
lap probability between exact and approximate states is 
0.997 in the worst possible case {n — 5, wq = i and 
ujfz = 25 Hz). In other words, the perturbation does not 
imply any significant deformation in the initial and final 
states. Its potential impact is instead at transient times 
where the energy may be high. 

We thus expand the overlap (?/'„(t/)|C/L(^/,0)|V'n(0)) 
as 1 + fn]n + /?vri + ... in powers of Vi to write 



Ft = 1+ f*^^ + f 



P(2) 



(A2) 



The first order correction is given by 



Jn,n ^ 



dt'{Mt')\Viit')\Mt')) 
'dt'L,l{t'){Mt')\AMt')), (A3) 



im 

and the unperturbed wavefunction corresponds to the 
known harmonic evolution of an expanding mode under 
the time-dependent frequency LUzjt) from an eigenstate 
of the initial harmonic oscillator 



e n 2b 

(62"n!)i/2 



Fin 



h b 



(A4) 

This is nothing but the eigenvector of the quadratic in- 
variant times the Lewis- Riesenfeld phase factor. 

Using the triangular inequality > ||x| — |?/||, 

Fl >l — \fn}n + fnX + ...|, and assuming that the per- 
turbative corrections satisfy \fn}n\ » |/ri^n| then 

FL>l-\fi]l\. (A5) 

The relevant matrix element in Eq. (jA3|) can be calcu- 
lated explicitly, 

2 



{Mt')\z'\Mt')) = 



3b^ 



^mujozj 4 
Using this result and the Ermakov equation, 



(A6) 



If (1)1 ^ 



[{r 



1) 



R 



1 



"^Qz Jo 



bb-^dt 



To bound — j^' bb^dt subjected to the imposed con- 
straints on 6 at t = and tf we follow closely the 
method applied in Q to bound the average energy. We 
integrate first by parts applying 5(0) = b{tf) — and 
minimize the resulting positive integral J^'^ b^b^dt, using 
the Euler-Lagrange equation bb + b^ ~ for &(0) = 1, 
b{tf) = 7. These two conditions are fulfilled for a set of 
functions larger than the subset of functions that satisfy 
all boundary conditions in Eq. (jl2p , so the minimization 
sets a lower bound for the later subset. The solution is 
b = [i(7^ — l)/tf + 1]^/^ and with it the integral becomes 



b^bMt = 



At 



f 



(A8) 



This b is not as such a physically valid solution, as the 
derivatives would jump at the edges, but it maximizes 
the fidelity lower bound. 



Fl > 1 



if 



3(7' - I f 



(A9) 
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Take note that there is an optimal t j value, and that the 
fidelity deteriorates when increasing the quantum num- 
ber and improves by increasing the waist. 

Using for b the quintic polynomial in Eq. (jl3p we get, 
instead of Eq. 



Now we apply this general expression to get the adi- 
abaticity condition on the longitudinal and radial fre- 
quencies of the Gaussian beam trap. In the harmonic 
approximation the longitudinal Hamiltonian is given by 
Eq. dini), 



h%^dt = ^"jl'_..^^' (iioi + 13517 + IIOI7'). 







24871t 



/ 



(AlO) 

This tends to OAA'y'^/tf for 7 >> 1, which, up to the 
constant factor, is the same dependence found for the 
optimal function. 

The first-order analysis provides analytical results and 
some qualitative guidance but for more accurate results 
we should resort to a second order calculation. Instead 
of Eq. (IA2p we may also write the fidelity for the n — th 
state in terms of the corresponding probability as Fl = 
(1 — X^n^^tn' Pn'Y^'^, whcrc P„' is the probability to find 
the system in the n'-th state a,t tf. This is approximated 
to second order using first order non-diagonal terms. 



Ft, 



\ 



p(i) 



(All) 



where 



/•(I) 

J n.n' 



dt'ilPnit'W^it'Mn'it')) 



(A12) 



/oo 
dy e-y"Hr,{y)H^,{y)y\ (A13) 
-00 



and 



Pn,n'{t) = / dh b^{ti)uji{ti)e 

Jo 



(A14) 



Appendix B: The condition of adiabaticity 

Consider a time-dependent Hamiltonian HQ{t), with 
instantaneous eigenstates and energies given by 

H,{t)m)^E,{t)m)- (Bi) 

The adiabaticity condition is given by 

\m\dtjm « ^mt)-E,{t)\, * ^ j. (b2) 



Hhar{z, t) — — 



2m dz^ 



■muj1{t)z'^ 



(B3) 



whereas the radial Hamiltonian is given, according to 
Eqs. (EH), (Eal), and by 



Hhar{r,t) 



2m dr^ 2 
2m 



(B4) 



The instantaneous eigenstates and energies of the Hamil- 
tonian, Eq. (|B3|) . are 



E„{t) = n = 0, 1, 



where uJz — i-^zit) and iJ„ is the Hcrmitc polynomial. To 
get the adiabaticity condition for the ground state i = 
of the longitudinal Hamiltonian we replace these expres- 
sions for the eigenstates and energies into Eq. (jB2[) . The 
overlap of the ground state with the first excited state 
J = 1 is so the first non- vanishing overlap of the ground 
state is with the j — 2 state. The adiabatic condition is 
satisfied if 



V2(^z 



< 1. 



(B5) 



For the radial Hamiltonian, Eq. (jB4[) . the instantaneous 
eigenstates and energies for 1/ = are (32| 



{r\Ht)) = 
Ek{t) = {2k 



2muJr 



-Lie- 

0,1,. 



k being the radial quantum number and the general- 
ized Laguerre polynomial. Again w/j = ujji{t). To get the 
adiabaticity condition for the radial direction we replace 



these expressions for j = and j = 1 into Eq. (|B2 



LOR 

4^1 



< 1. 



(B6) 
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